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We report on the status of our study towards the equation of state in 2 + 1 flavor QCD with 
improved Wilson quarks. To reduce the computational cost which is quite demanding for Wilson- 
type quarks, we adopt the fixed scale approach, i.e. the temperature T is varied by A^, at fixed 
lattice spacing. Since the conventional integral method to obtain the pressure is inapplicable 
at a fixed scale, we adopt the "T-integral method", to calculate the pressure non-perturbatively. 
Reduction of the computational cost of T = simulations thus achieved is indispensable to study 
EOS in QCD with dynamical quarks. 
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1. Introduction 

Clarification of tiie equation of state (EOS) of iiot QCD is important to understand the nature 
of quark matter in early Universe and in relativistic heavy ion collisions. Most lattice studies have 
been done with computationally cheap staggered-type lattice quarks. However, their theoretical 
basis such as locality and universality are not well established. Therefore, to evaluate the effects of 
lattice artifacts, it is important to compare the results with those obtained using theoretically sound 
lattice quarks, such as the Wilson-type quarks. 

In this note, we report on the status of our study towards the EOS in QCD with 2 + 1 flavors 
of dynamical Wilson-type quarks. To reduce the lattice artifacts, we adopt RG-improved Iwasaki 
gauge action and the clover-improved Wilson quark action with non-perturbatively adjusted clover 
coefficient. 

A reason that Wilson-type quarks have not been intensively studied in finite temperature QCD 
is that the computational cost for Wilson-type quarks is larger than that for staggered-type quarks, 
in particular at small quark masses. Therefore, we have to implement efficient methods for sim- 
ulations and analyses. We adopt a fixed scale approach in which the pressure is calculated non- 
perturbatively by the T-integral method 

2. Fixed scale approach armed with the T-integral method 

Conventionally, finite temperature simulations in lattice QCD are performed in the fixed-A/, 
approach, where temperature T = {Nfa)^^ is varied by changing the lattice scale a (through a 
variation of the lattice gauge coupling j8) at a fixed temporal lattice size A^,. Thus, simulations 
have to be repeated at different values of j8 along a line of constant physics (LCP) in the coupling 
parameter space. In this approach, a sizable fraction of the computational cost is devoted for T = 
simulations to set the basic parameters such as the lattice scale, to determine LCP's and the beta 
functions on them, and to carry out zero-temperature subtractions for the renormalization of finite- 
temperature observables. Note that these zero temperature simulations are required at all the points 
in the coupling parameter space for finite temperature simulations. 

In the fixed scale approach we adopt, we vary T by changing Nt at a fixed a, fixing all coupling 
parameters. Since the coupling parameters are common to all temperatures, (i) T = subtractions 
can be done by a common zero temperature simulation, (ii) the condition to follow the LCP is 
obviously satisfied, and (iii) the lattice scale etc. are required only at the simulation point. We may 
even borrow results of existing high precision spectrum studies at T = which are public e.g. on 
the International Lattice Data Grid (ILDG) Then, the computational cost needed for T = 
simulations can be reduced largely. 

Because the lattice spacings in spectrum studies are usually smaller than those used in con- 
ventional fixed-A/f studies around the critical temperature T^, for thermodynamic quantities around 
Tc, we can largely reduce the lattice artifacts due to large a and/or small Nt than those in the con- 
ventional fixed-A/, approach. This requires a larger computational cost at low temperatures due 
to the larger A',. Nevertheless, the merits around Tc will be a good news for phenomenological 
applications, since temperatures relevant at RHIC and LHC will be at most up to a few times Tc. 
On the other hand, as T increases, A'f in our approach becomes small and hence the lattice artifact 
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increases. Therefore, our approach is not suitable for studying the high temperature Umit. Note 
that the pros and cons of our method are complement to the conventional method. 

In the conventional fixed-A/, approach, p is calculated non-perturbatively by the integral method 
[^]: using the thermodynamic relation p = {T /V) InZ valid in the large volume limit, with V being 
the spatial volume and Z the partition function, 

where S is the action and b = {p, K^d, Ki, • " ■) = {bijbx,- ■■) is the vector in the coupling parameter 
space. The integration path can be chosen freely in the coupling parameter space as far as p(^o) ~ 0- 
The conventional integral method is inapplicable in the fixed scale approach because we sim- 
ulate only at a point in the coupling parameter space. Therefore, we developed "the T-integral 
method" to evaluate the pressure non-perturbatively: Using a thermodynamic relation vahd at 
vanishing chemical potential 

±(J^\e-3p 

dT\T^) r4 ' ^ ' 



we obtain 



with p{To) ~ 0. Here the trace anomaly e — 3pis calculated as usual by 

3p _ y dbi f / dS\ / dS\ 



-N^l^''ld[\dbJ-\dbXJ ^'-'^ 

where A^^ is the spatial lattice size. The coefficient is the lattice beta-function defined by the 
variation of the ith coupling parameter bi along the LCP. 

In the fixed scale approach, T is restricted to discrete values due to the discreteness of A^,. For 
the integration of (^!3|), we need to interpolate the data with respect to T. The systematic error from 
the interpolation should be checked. Note that, because the scale is common for all data points in 
the fixed scale approach, T is determined without errors besides the common overall factor I /a. 



3. Test of the method in quenched QCD 

In |[T]], we made a test of the fixed scale approach armed with the T-integral method in quenched 
QCD using the standard one plaquette gauge action, on isotropic and anisotropic lattices. The 
simulation parameters are summarized in Table |I[ 

The trace anomaly obtained on the isotropic il, 12 and 13 lattices are shown in the left panel 
of Fig.|l|. The shaded line represents the result of the conventional fixed-A^, method obtained on a 
large lattice of N, = S and A^, = 32 (about 2.7 fm around Tc ~ 290 MeV) [|]. Comparing il and 
13, we find that the lattice cutoff effects are quite small on these lattices. On the other hand, the 12 
lattice shows a small deviation from the Nt = S lattice near Tc. This deviation may be understood 
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set 


iS 






N, 






L[fm] 


a{dg '^/da) 


il 


6.0 


1 


16 


3-10 


5.35(^2) 


0.093 


1.5 


-0.098172 


i2 


6.0 


1 


24 


3-10 


5.35(^2) 


0.093 


2.2 


-0.098172 


i3 


6.2 


1 


22 


4-13 


7.37(3) 


0.068 


1.5 


-0.112127 


a2 


6.1 


4 


20 


8-34 


5.140(32) 


0.097 


1.9 


-0.10704 



Table 1: Simulation parameters on isotropic and anisotropic lattices The il, i2 and i3 lattices are 
isotropic, while the a3 lattice is anisotropic with ^ = as /at = 4. The beta function is taken from [^. The 
temperature ranges cover T ^ 200-700 MeV. Corresponding T = simulations are done on A^, =20*^ 
lattices. 
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Figure 1: EOS in quenched QCD [jl[|. Left: trace anomaly on isotropic lattices. Dotted lines are the natural 
cubic spline interpolations of the data. Center: trace anomaly on anisotropic lattice a2 compared with the 
isotropic i2 lattice with similar spatial lattice spacing and volume. Right: energy density and pressure by 
the T-integral method. The shaded curves represent the results of the conventional fixed-A^( method at A^, = 8 




by the physical finite size effect expected in the critical region. Off the critical region, all results 
agree well with each other. 

Dotted curves in the left panel of Fig.|] are the natural cubic spline interpolations of our trace 
anomaly. To estimate the systematic error due to the interpolation ansatz, we tested another inter- 
polation with the trapezoidal rule. Carrying out the numerical integration (2.3), we find that the 
EOS from the trapezoidal interpolation is almost identical with the EOS from the natural cubic 
spline interpolation [[T]]. 

To further estimate systematic effects due to discreteness of T, we compare the results with 
those on the anisotropic lattice a2, which has about 4 times finer resolution in T than the 12 lattice. 
In the central panel of Fig.[l|, we compare the trace anomaly on a2 and 12 lattices. We find that the 
data points from the a2 lattice are well on the natural cubic spline interpolation line of the 12 lattice, 
except for the data on the a2 lattice around the peak where the interpolation line of the 12 lattice 
slightly overshoots. We note that the height of the peak on the a2 lattice is consistent to those of 
the fine 13 and Nt = S lattices shown in the left panel of Fig.|l]. Therefore, the difference may be 
explained by the smaller discretization errors in the temporal direction on the a2 lattice. 

In the right panel of Figjl], we show the results of EOS by the numerical integration (2.3). We 
find that the results of 12 and a2 lattices are well consistent with each other. This suggests that the 
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Figure 2: Status of finite temperature 2+1 flavor QCD simulations with improved Wilson quarks. Left: 
the Polyakov loop time history. Right: heavy quark free energy in the color-singlet channel [pO- The 
heavy-quark potential V{r) at T = was calculated by the CP-PACS and JLQCD Collaborations [^]. The 
pale dashed curve is the result of a Coulomb + linear fit of y(r). The arrows on the right side denote twice 
the thermal average of the single-quark free energy. 



systematic errors due to the discreteness of T is at most about the statistical errors. The shaded 
curve in the figure represents the result of the conventional fixed-Nt method at Nt = S [Q]. We 
find that the fixed scale approach armed with the T-integral method is powerful to calculate EOS 
reliably. See ReU^ for more discussions. 

4. Towards the EOS in 2 + 1 flavor QCD with improved Wilson quarks 

Adopting the fixed scale approach armed with the T-integral method, we are carrying out a 
series of simulations of finite temperature QCD with 2 + 1 flavors of improved Wilson quarks. 
As the basic zero temperature configurations, we use those created by the CP-PACS and JLQCD 
Collaborations [^] and made public at the JLDG branch of ILDG [Q] . Their spatial lattice volume 
is about (2 fm)"*. Among their simulation points, we have chosen the finest lattice (a = 0.07 fm, 
j8 = 2.05) with the lightest u and d quarks {rrij^/mp = 0.63) The lattice size is 28^ x 56 and the 
statistics is about 6000 trajectories. Using the same coupling parameters as the zero temperature 
simulation, we are generating finite temperature configurations on 32^ x Nt lattices with A^, = 4, 6, 
• ■ •, 16. The pseudo-critical temperature is expected to be around Nt ~ 14. 

Status of the finite temperature simulations is shown in the left panel of Fig.^. While we 
are still increasing the statictics, in particular for Ni = 12-16 lattices around the pseudo-critical 
temperature, we have started first test calculations on these configurations. 

At the conference, Yu Maezawa presented our preliminary results for the heavy quark free 
energy obtained on these configurations [^. Our results for the heavy quark free energy in the color 
singlet channel are shown in the right panel of Fig.^ Here, we note another good feature of the 
fixed scale approach that we can purely extract the temperature effects in the physical observables: 
In the case of fixed-A/; approach, because the j3 -dependent renormalization factor for the constant 
term of the free energy is not known, the free energies at different temperatures (different j8 's) are 
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Figure 3: mps/my ratio in 2 + 1 flavor QCD. Left: mps/my for light-light i£ mesons (m,i/mp). Right: 
mps / my for ss mesons. Curves represent the results of a 9 parameter global fit. 



vertically adjusted by hand such that they coinside with each other at short distances. This means 
that we imply the temperature dependence to be small at short distances. On the other hand, in 
the fixed scale approach, the renormalization factors are common to all temperatures. Therefore, 
we need no adjustment of the constant term to compare the free energies at different temperatures. 
The free energies shown in the right panel of Fig.^are plotted without adjusting the constant term, 
and V{r) is the zero-temperature heavy quark potential defined by Wilson loop expectation values 
We find that the free energies at various T converge to V{r) at short distances. We have thus 
proved the validity of the theoretical expectation that the short distance physics is insensitive to the 
temperature. See [H] for more discussions. 



We now turn our attention to the calculation of EOS. According to ( |2.4[ ), in addition to the 
gluon contribution to the trace anomaly, we have the quark contribution due to the scaling of the 
hopping parameters K"„j and K^. When the clover coefficient Csw depends on j3, its j3 -derivative 
also contributes as a part of the quark contribution. Therefore, we need to know the values of the 
beta functions for these coupling parameters. 

We attempt to calculate the beta functions in 2 + 1 flavor QCD using the hadron data by the 
CP-PACS and JLQCD Collaborations In this study, we use the data of mps/niy for light-hght il 
mesons {nij^/mp), mps/my for ^'5 mesons, and the decay constant /ps of the ^'^^pseudoscalar meson, 
to obtain the LOC through our simulation point mps/mv(^^) = 0.6337 and mps/my{ss) = 0.7377 
in the three dimensional coupling parameter space of (jS, K^d, Ks)^ as well as the scale on the LCP 
Figure || shows the data of mps/my{£i) and mps/mv(i^). The curves in the plots are the results of 
a 9-parameter global fit. Although the fit approximately reproduces the data, the quality of the fit 
is not quite high with Mof ~ 10. To calculate the beta functions, we adopt the inverse matrix 
method developed in [^. Results of the beta functions are shown in the left panel of Fig.^. We find 
that the magnitudes of the beta functions are similar to those obtained in the previous two flavor 
case with improved Wilson quarks [^. However, although the beta function adp /da for the gauge 
coupling is well determined, the beta functions adKud/da and adKg/da for the hopping parameters 
have errors larger than their values. The main reason is the coarseness of the data points in the 
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Figure 4: Beta functions and EOS in 2 + 1 flavor QCD. Left: beta functions adji/da, adKi,d/da and 
adKs/da. Right: preliminary result of the gluon contribution to the trace anomaly. 



coupling parameter space. The substitution of other meson masses for /ps (ss) leads to similar 
results. 

Using the result of a dji /da, we further attempt to calculate the gluon contribution to the trace 
anomaly. A preliminary result is shown in the right panel of Fig.^. For comparison, the peak height 
of {e — 3p)/T'^ in the case of 2 -|- 1 flavor QCD with improved staggered quarks was about 6-8 on 
Nt = 6 and 8 lattices On the other hand, we expect a large cancellation between the gluon and 
quark contributions in the trace anomaly: In the case of two flavor QCD with a similar improved 
Wilson quark action ^, the peak height of {e — 3p)/T'^ is about 13 at mji/nip ~ 0.65 on A'f = 4 
lattices, in which the gauge contribution is about 45 and the quark contribution is about —32. Thus 
the magnitude of the gluon contribution shown in the right panel of Fig.^ seems to be consistent 
with expectation. 



5. Discussion 

We have developed the fixed scale approach armed with the T-integral method to reduce the 
computational cost for the EOS calculation on the lattice. A test in quenched QCD has shown that 
the method works well 

Applying the method, we are carrying out a series of finite temperature simulations in 2 -|- 1 
flavor QCD with improved Wilson quarks, based on the pubUc zero-temperature configurations on 
ILDG generated by the CP-PACS and JLQCD Collaborations [g]. 

To calculate the EOS, we need the beta functions too. We found that the hadron data available 
from the zero-temperature spectrum study are insufficient to calculate precise beta functions for the 
hopping parameters — we need more data points in the coupling parameter space around the finite 
temperature simulation point. In order to avoid additional intensive zero-temperature simulations, 
we are trying to develope a reweighting technique to directly calculate the beta functions at the 
simulation point. 
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